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ABSTRACT 

The  principal  problem  associated  with  steady-state  simulation  is  the 
estimation  of  the  variance  term  in  an  associated  central  limit  theorem. 
This  paper  develops  several  strongly  consistent  estimates  for  this  terra 
using  the  strong  approximations  available  for  Brownian  motion.  A 
comparison  of  rates  of  convergence  is  given  for  a  variety  of  estimators. 
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1 


INTRODUCTION 


Let  X  *  {X(t)  :  t  >_  0}  be  a  real-valued  (measurable)  stochastic  pro¬ 
cess  representing  the  output  of  a  simulation.  (To  incorporate  output 

sequences  X  into  our  framework,  we  set  X(t)  *  X.  ,  ,  where  (t|  is  the 
n  1 1 1 

greatest  integer  less  than  or  equal  to  t.)  The  process  X  is  said  to 
possess  a  steady-state  if  there  exists  a  finite  constant  r  such  that 

,  t 

(1.1)  r(t)  =  7  J  X(s)ds  — »  r 

c  0 

as  t  ♦  where  denotes  weak  convergence.  The  problem  of  consistent¬ 

ly  estimating  and  producing  confidence  intervals  for  the  parameter  r  is 
known,  in  the  simulation  literature,  as  the  steady-state  simulation  prob¬ 
lem. 

The  limit  theorem  (1.1)  suggests  that  r(t)  can  be  used  to  consis¬ 
tently  estimate  r.  It  turns  out  that  one  can  frequently  establish  that 
the  output  process  X  in  fact  satisfies  a  somewhat  stronger  relation, 
namely,  there  exist  finite  constants  r  and  a  such  that 

(1.2)  t1/2(r(t)  -  r)  ^  a  N(0,1)  . 

Suppose  now  that  one  can  construct  an  estimator  s(t)  such  that 

(1.3)  s( t)  o 

as  t  +  ®.  Then,  (1.2)  and  (1.3)  together  imply  that  if  o  >  0,  then  the 
interval 


l 


(1.4) 


[r(t)  -  z(6) 


s(t) 


r(t)  +  z( 6) 


is  an  asymptotic  100(1-6)%  confidence  interval  for  r,  provided  that  2(6) 
is  chosen  to  solve  the  equation  P(N(0,l)  _<  z(  6) }  ■  1  -  6/2. 

The  above  discussion  suggests  that  in  the  presence  of  an  output  pro¬ 
cess  satisfying  (1.2)  with  a  >  0,  the  steady-state  simulation  Droblem  is, 
in  principle,  solved,  once  an  estimator  s(t)  for  a  has  been  constructed. 
Thus,  the  determination  of  an  estimator  s(t)  for  a  can  be  viewed  as  the 
fundamental  problem  of  steady-state  simulation  output  analysis.  (It 
should,  however,  be  noted  that  certain  output  analysis  methods  employ  a 


s 


M 


different  approach,  which  does  not  require  explicit  estimation  of  0;  the 
methods  of  batch  means  (see  pp.  85-89  of  BRATLEY,  FOX,  and  SCHRAGE  (1983) 
for  a  description)  and  standardized  time  series  (SCHRUBEN  (1983))  fall  into 
this  category). 

As  a  consequence  of  the  above  observation,  a  great  deal  of  attention 
has  been  focused  on  the  construction  of  such  estimators  for  o.  In  some 
sense,  all  currently  available  estimation  methods  make  use  of  the  fact  that 
if  X  is  well-behaved  and  approximately  stationary,  then  (1.2)  suggests 
that 


lim  E(t(r(t)  -  r))2  -  o2EN(0,l)2 


(1.5) 


i.e.,  2  /  EX  ( 0)x  (t)dt  -  a2 

0  C  c 


M 


i 


a 


where  Xc<t)  =  X(t)  -  r.  Spectral  procedures  (see  pp.  95-98  of  [3D  use 
kernel  estimators  to  estimate  the  left-hand  side  of  (1.5),  whereas  autore¬ 
gressive  methods  (see  pp.  98-101  of  { 3 J )  fit  an  autoregressive  process  to 
X,  compute  the  analog  of  (1.5)  for  the  fitted  process,  and  use  that  quan¬ 
tity  as  an  estimator  of  (1.5)  for  X.  As  for  the  regenerative  method  (see 
pp.  89-94  of  [3D,  it,  loosely  speaking,  uses  the  special  structure  of  a 
regenerative  process  to  "truncate”  the  integral  on  the  left-hand  side  of 
(1.5)  at  a  regeneration  time,  threby  simplifying  the  estimation  problem. 

In  this  paper,  we  propose  a  new  method  for  consistently  estimating  a, 
which  does  not  make  any  explicit  use  of  relationship  (1.5).  Our  estimators 
are  based  on  known  limit  theorems  about  the  increments  of  Brownian  motion; 
the  constant  a  appears  as  a  scaling  constant  in  these  limit  theorems.  We 
then  use  strong  approximation  methods  to  translate  the  resulting  estimator 
for  a  from  the  Brownian  motion  back  to  the  original  output  process  X; 
the  resulting  estimator  for  a  depends  only  on  the  observed  values  of  X, 
and  not  the  Brownian  motion. 

In  some  respects,  our  method  is  similar  to  that  of  SCHRUBEN  (1983). 

The  method  of  standarized  time  series  depends  on  the  fact  that  a  appears 
as  a  scaling  constant  when  one  weakly  approximates  the  original  process  by 
a  Brownian  motion;  a  is  not  estimated  but  is  instead  "cancelled"  out.  Bv 
contrast  our  method  involves  strong  approximation  results  and  gives  rise  to 
strongly  consistent  estimators  for  a.  As  proved  in  GLYNN  and  IGLEHART 
(1985),  consistent  estimation  of  o  enjoys  certain  asvmptotic  advantages 
over  standardized  time  series. 

In  Section  2,  we  introduce  our  estimators  s(t)  for  a;  our  basic 
hypothesis  on  X  is  that  a  suitable  strong  approximation  by  Brownian 


motion  is  possible.  Section  3  is  devoted  to  discussion  of  processes  X 
which  satisfy  the  strong  approximation  hypothesis.  In  Section  4,  the  rate 
of  convergence  of  s(t)  to  o  is  studied,  and  compared  to  that  available 
via  the  regenerative  method.  Our  main  contribution  here  is  to  suggest  that 
entirely  new  methods  for  estimating  a  may  be  worth  exploring.  Further 
comparison  of  these  new  methods  with  the  old  methods  should  be  carried  out 
via  numerical  examples. 

2.  A  NEW  CLASS  OP  ESTIMATORS  FOR  <7 

Let  S  *  (S(t):  t  _>  0) ,  where  S(t)  *  ^  X(s)ds.  Throughout  this 
section,  we  shall  assume  that: 

(2.1)  there  exists  a  probability  space  supporting  a  process  S*  and  a 

standard  Brownian  motion  such  that: 

(i)  the  distribution  of  S  equals  that  of  S* 

1  /?— X 

(ii)  S*(t)  -  rt  »  oB(t)  +  0(t  )  a.s.  for  some  constants  r, 

o,  and  X  (0  <  X  <  1/2,  a  >  0)  as  t  ♦  •. 

We  shall  refer  to  (2.1)  as  our  strong  approximation  assumption;  It  says 
that  a  process  S*,  possessing  precisely  the  same  distribution  as  S,  can 
be  a.s.  well  approximated  by  a  Brownian  motion.  Mote  that  under  (2.1) 

(ii), 

(2.2)  t_1/2(S*(t)  -  rt)  -  t"1/2o  B(t)  ♦  0  a.s. 

as  t  ♦  ••  By  (2.1)  (i),  it  follows  that  (2.2)  also  holds  with  S  taking 
the  role  of  S*;  this  shows  that  (2.1)  automatically  implies  (1.2). 


For  0  <  p  <  1,  let 

s  j  ( t)  -  sup  S(u  +  tP)  -  r(t)tP  ~  S(u) 

0<u<t-tP  ( 2tP  *  (1  -  p)  *  log  t)1/2 

(2.3)  THEOREM.  If  p  Is  chosen  so  that  p+2\>  1,  0<p<  1,  then 
s^(t)  ♦  or  a.s.  as  t  ♦  «. 

PROOF.  CSORGO  and  REVESZ  (1979a)  showed  that 

(2.4)  lim  sup  B(u  +  tP)  -  B(u)  ,  _ 

t+"  °<u£c“cP  (2tp  *  (l  -  p)  *  io«  t)1/2 

Let  S*(t)  -  S*( t)  -  rt.  By  (2.1)  (ii),  S*(t)  *  OB(t)  +  (Xt1^2  S  a.s. 
c  c 

so  it  follows  that 


(2.5) 


sup  |oB(u  +  tP)  -  <JB(u)  -  S*(u  +  tP)  -  S*(u)  ! 

c  c 

0<u<t-tP 


-  0(t 


1/2-^ 


a.s. 


Relations  (2.4)  and  (2.5),  together  with  the  condition  2p  +  \  >  1,  imply 
that 


11m 

sup 

S*(u  +  tp)  -  S*(u)  a 

c  c  *  o  a.s. 

t-*“ 

0<u<t-tP 

( 2 1 p  •  (1  -  P)  log  t)1/2 

lim 

sup 

S*(u  +  tp)  -  rtp  -  S*(u)  ■  a  a.s. 

CKu<t-tP 

(2tp  •  (1  -  p)  log  t)1^2 

Furthermore, 

the  law 

of  the  iterated  logarithm  for  Brownian  motion 

implies  that 


Applying  the  strong  approximation  (2.1)  (ii)  to  (2.7),  we  find  that 


(2.8) 


r*(t)  -  r 


0  ^ (log,  log  tjl/2) 


3  •  S  • 


where  r*(t)  *  S*(t)/t.  Since  (log  log  t)1/2  •  t(p-1)  >  0,  it  follows  from 
(2.6)  and  (2.8)  that 


(2.9) 


lim  sup 
t+"  0<u<t-tp 


S*(u+tP)  -  tPr*( t)  ~  S*(u) 
(2tP  •  (1-p)  •  log  t)1/2 


@  3  •  S  • 


But  S  has  the  same  distribution  as  S*,  so  the  theorem  follows  immedi¬ 
ately  from  (2.9). 


We  can  further  define  the  following  estimators: 


s2(t)  -  sup 

0<u<t-tP 


S(u+tP)  -  r(t)tP  -  S(u) 

( 2 1 P  •  (1-p)  *  log  t)1/2 


a(t)  -  (8  •  (1-p)  •  log  t/n2cp)1/2 

(2.10)  THEOREM.  If  p  is  chosen  so  that  p+2X>  1,  0<p<  1,  then 
s^t)  +  a  a.s.  as  t  ♦  ®,  2  <  i  <  5. 

The  proof  of  this  theorem  follows  exactly  as  for  Theorem  2.4;  the  key 
is  to  have  available  an  analog  to  (2.4)  for  Brownian  motion,  for  each  of 
the  estimators  s^(t),  2  _<  i  £  5.  For  s^(t),  2  £  i  £  4,  we  refer  to  CSORGO 
and  REVESZ  (1979a);  for  s^(t),  the  result  can  be  found  in  CSORGO  and 
REVESZ  (1979b). 

3.  STRONG  APPROXIMATION  OP  STOCHASTIC  PROCESSES 

We  now  proceed  to  discuss  conditions  under  which  Assumption  2.1  holds. 

SETTING  Is  Let  T  ■  (Y(t)  :  t  ^  0)  be  a  (possibly)  delayed  regenerative 
process  with  regeneration  times  0£  T(0)  <T( !)<•••  If  f  is  a  real- 
valued  function  defined  on  the  state  space  of  T,  set  X(t)  *  f(Y(t)).  Let 
x(k)  -  T(k)  -  T(k-l)  for  k^  l,  and  assume  that  Y  is  positive  recurrent 
in  the  sense  that  Et(l)  <  «.  Suppose  further  that  there  exists  0  <  6  <  2 


such  that 


(X(s)  -  n)  ds  and 


(ill)  EZ2(1)  >  0,  where  Z(n)  -  / 

T(n-l) 


T(  1 ) 

H  =  E (/  X(s)  ds)/Ex(l)  . 

T(0) 


Then,  (2.1)  is  satisfied  with  r  -  n,  a2  -  EZ2( 1 ) /Ex(l ) ,  and  X.  satis¬ 
fying  0  <  \  <  6/(4  +  26);  see  pp.  117-122  of  PHILIPP  and  STOUT  (1975)  for 
a  proof  in  the  case  where  T  is  a  countable  state  irreducible  Markov  chain 
(their  argument  easily  adapts  to  the  more  general  regenerative  setting 
described  above). 


SETTWG  2:  Let  X  •  {x^  :  n  _>  0)  be  a  strictly  stationary  sequence  of 

a  i  I 

r.v.'s  for  which  there  exists  0  <  6  <  1  such  that  E | XQ |  <  «.  Sup¬ 

pose,  in  addition,  that  X  is  ♦-mixing  (see  Section  20  of  BILLINGSLEY 
(1968)  for  a  discussion)  with  mixing  coefficients  satisfying 
,  *(n)1/2  <  If 

n- 1 


(3.2) 


a  -  2  /  EX  (0)X  (s)ds  >  0 

A  c  C 


(the  integral  (3.2)  converges  absolutely),  then  (2.1)  is  satisfied  with 
r  ■  EX( 0) ,  o2  •  a,  and  X  satisfying  0  <  \  <  6/(24  +  126)  (if  XQ  is 
bounded  a.s.,  then  X.  can  be  chosen  to  be  1/12).  This  result  can  be 
found  on  pp.  26-38  of  [11].  (For  a  version  of  this  result  in  the  case  when 
I*  ,  $(n)1/2  <  -,  see  BERKES  and  PHILIPP  (1979).) 


Farther  strong  approximation  theorems  are  also  available  for  lacunary 
trigonometric  series,  martingale  processes,  Gaussian  sequences , and  strongly 


mixing  processes;  see  (111  for  a  complete  description  of  such  results. 

Thus,  the  assumption  (2.1)  is  satisfied  by  a  large  class  of  stochastic  pro¬ 
cesses  exhibiting  weak  time  dependencies. 

We  should  also  comment  that  the  convergence  rates  (i.e.,  the  size  of 
X)  quoted  above  for  regenerative  processes  and  ^-mixing  sequences  can 
probably  be  improved.  For  example,  much  better  results  are  available  for 
sequences  of  independent  and  identically  distributed  (i.i.d)  r.v.'s.  In 
particular,  it  is  reasonable  to  expect  that  results  for  regenerative 
processes  can  be  obtained  for  X  arbitrarily  close  to  1/2. 

SETTING  3:  Let  X  -  (Xn  :  n  >  0>  be  a  sequence  of  non-degenerate 
i.i.d. r.v.’s  with  E|Xq|p  <  «  for  p  >  2.  Then,  (2.1)  is  satisfied  with 
r  -  EXq,  o2  -  var(Xp) ,  and  X  *  1  -  1/p;  see  KOMLOS,  MAJOR,  and  TUSNADY 
(1975,1976)  and  MAJOR  (1976)  for  proofs. 


4.  COMPARISON  OF  CONVERGENCE  RATES 

It  has  been  shown  by  REVESZ  (1980)  that  if  0  <  p  <  1,  then 


lim  (log  t) 


i*6 


B(u+tP)  -  B(u) 

sup  - - - pry 

(Xu<t-tP  (2t  *  (l'p)  ’  l0*  C) 


-  a 


*  0  a .  s . 


lim  (log  t) 

t+m 


7  +6 


sup  sup 


B ( u+v )  -  B(u) 


0<u<t-tP  (Xv<tp  (2tP  *  (l_p)  *  loa  C) 


777 


-  a 


0  a.s. 


for  6  -  0;  it  is  further  Indicated  in  CSORGO  and  STEINEBACH  (1981)  that 
for  6  >  0,  the  above  lim  sup's  are  infinite. 


Using  the  strong  approximation  (2.1),  it  follows  that 


lim  (log  t)d/2)+6  s  (t)  -  a  -0  a.s. 

t-*.® 


for  6  *  0  (i  *  1,3),  whereas  divergence  occurs  if  6  >  0.  Thus,  the  rate 


convergence  of  s^t)  (i  «  1,3)  to  a  is,  roughly  speaking,  of  order 


Uo, 


It  is  instructive  to  compare  this  rate  of  convergence  to  that  avail¬ 


able  when  a  is  estimated  via  the  regenerative  method  of  simulation  (we 


choose  this  method  as  a  basis  for  comparison,  since  we  can  do  the  conver¬ 


gence  rate  analysis  easily  in  this  setting). 


Let  T  be  a  regenerative  process  with  regeneration  times  0  <  T(0) 


<  T(l)  <  •••;  set  X(t)  *  f(Y(t)),  where  f  is  a  real-valued  function 


defined  on  the  state  space  of  T.  Put  T(-l)  *  0  and  let  N(t)  - 


max  (k  >  -1  :  T(k)  <_  t).  The  basic  regenerative  estimator  for  a  is  given 


0  ;  N(t)  <  0 


[{  1  (V1  -  r(t)  x^2]172  ;  N(t)  >  1  , 


where 


V  -  /  X(s)ds  and  x  -  T(i)  -  T(i-l). 

T(i-l)  1 


'  vv->vywu 


as 


(4.2)  THEOREM.  If  E (/  ( !x(s)  |  +  l)ds)  <  •, 

T(i-l) 

|s(t)  -  o  |  -  a.s. , 

where 

3  -  — _i -  •  e[A.  -  Xz.)2  , 

•  rt 

Z  -  V  -  rx  , 
n  n  n 

2  2 

A  *  Z  -at  ,  and 
n  n  n 

X  -  2  EZjTj/ETj  . 

Recall  that  for  regenerative  processes,  r  ■  EV^/Ex^  and  o2  * 

2 

EZ./Ex  ,  so  that  Z  and  A  are  mean-zero  r.v.'s;  we  will  need  this 
11  n  n 

fact  in  our  proof  of  Theorem  4.2.  Also,  we  remark  that  Theorem  4.2  is  a 
statement  of  the  law  of  the  iterated  logarithm  for  the  estimator  s(t). 

PROOF  (of  Theorem  4.2):  On  the  event  (N(t)  l),  observe  that  if  v(t) 

2 

s  (t) ,  then 


m 


§ 


(4.3) 


2  1  N£C)  2 
v(t)  -  a2  -  i  J  Z2 

i-1  1 


1  N£c) 

+  2(r  -  r(  t) )  •  ±  I  Z1xi 


+  (r  -  r(t))‘ 


•  i  1°  x2 
t  A.  i 


.  N(t) 

t  l  (Ai  - 

C  i-1  1  1 


t  N(t) 

-2  /  (X(s)  -  r)ds  •  I  Z  x 

T(N(t) )  A  i-1 


.  N(t)  EZ.T 

2  (r  -  r(c))  •  (i  j  Vl’-iT1) 


2  1  2 
(r  -  r(t))2  •  ±  2  xf  . 

t  »_  i  1 


But  s(t)  ■  g(v(t))  where  g(x)  *  x1  f  so  by  Taylor* s  theorem,  we  have 


s(t)  -  a  +  g‘(v(t))  (v(t)  -  o2) 


,  2  l  /? 

where  5(t)  lies  between  v(t)  and  a  and  g’(x)  -  l/(2x  ).  Since  v(t) 


o  a.s.,  it  follows  that  g’(5(t))  ♦  l/(2o)  a.s.  Thus,  to  prove  the 


theorem,  it  suffices  to  show  that 


(4.4) 


IS  V  2~Iog^Tog  t  I  v(t)  -  g2  j  “  2°P1/2  a.s. 


V 

I 

s 


B 


By  the  Hartman-Wintner  law  of  the  iterated  logarithm  and  the  fact  that 


a 


But  N(t)/t  +  1/E't1  a.s.  as  t  ♦ 


(4‘5)  V  TToglog 


►  I  .  N(c) 

- -  -  l 

log  t  it  1_1 


(4-6)  iog-io8 1  <r  - r<t)) 


(4.7)  lim 


log  log  t 


(r  -  r(t) )' 


Furthermore, 


;n 


t  T(k) 

/  |X(s)  -  r|  ds  _<  max  (/  |x(s)  Ids  +  |r|  t  )  . 

T(N(t) )  l<k<N(t)+l  T(k-l) 


Our  moment  hypothesis  allows  one  to  apply  the  Borel-Cantelli  lemma  to 


obtain 


as  k  ♦  •;  this  shows  that 


T(k) 

max  (/  !x(s) Ids  +  |r |  t  )  >  0  a.s . 

l<k<N(t)+l  T(k-l)  K 


and  t  ♦  ®.  Hence 


(4.8) 


_ _  t  N(t) 

lim  tJ/  /  (X(s)  -  r)ds  •  I  Z  x 

t-*»  T(N(  t) )  t  i-1 


0  a.s . 


as  f*®.  Combining  (4.5),  (4.6),  (4.7),  and  (4.8),  we  see  that  the  decom¬ 
position  (4.3)  yields  (4.4). 

Roughly  speaking,  Theorem  4.2  says  that  s(t)  converges  to  o  at 
1/2 

rate  (log  log  t/t)  .  By  comparison  with  the  previously  obtained  con¬ 
vergence  rate  for  3^(0  (i»l,3),  this  is  much  faster.  This  does  not 
necessarily,  however,  imply  that  the  estimators  3^(0  for  o  will  behave 
worse  than  the  regenerative  estimator  for  purposes  of  confidence  interval 


generation 
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